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MULTIPLE SHAPE REGISTRATION USING CONSTRAINED OPTIMAL CONTROL 


SYLVAIN ARGUILLERE, EMMANUEL TRELAT, ALAIN TROUVE, AND LAURENT YOUNES 


Abstract. Lagrangian particle formulations of the large deformation diffeomorphic metric mapping algorithm (LDDMM) 
only allow for the study of a single shape. In this paper, we introduce and discuss both a theoretical and practical setting for 
the simultaneous study of multiple shapes that are either stitched to one another or slide along a submanifold. The method 
is described within the optimal control formalism, and optimality conditions are given, together with the equations that are 
needed to implement augmented Lagrangian methods. Experimental results are provided for stitched and sliding surfaces. 


1. Introduction 

The large deformation diffeomorphic metric mapping (LDDMM) approach to shape matching is a powerful topology¬ 
preserving registration method with an increasing record of successful applications in medical imaging. It was 
first described in 1381 for point sets and in im EH Hi la for images and has become widely used in the med¬ 
ical imaging literature and other applications. While deeper understanding and extensions of the underlying the¬ 
oretical framework was pursued El Ea uni Ea cni ED El iH la and alternative numerical methods were designed 
|5l[T4l|60l|43|22llMll31[23[I3 LDDMM has been applied to medical imaging data including brain 14511771[5^fTTlIS?! . 
heart mis and lung l65l images. This algorithm provides a non-rigid registration method between various types of 
objects (point sets, curves surfaces, functions, vector fields...) within a unified framework driven by Grenander’s 
concept of deformable templates 1^ . It optimizes a flow of diffeomorphisms that transform an initial object (shape) 
into a target one. 

The practical importance of shape registration is underlined by the increasing amount of work that has flourished 
in the literature over the past few years. LDDMM is one among many methods that have been proposed to perform 
this task. Several such methods are based on elastic matching energies Qna, and other, like LDDMM, inspired by 
viscous fluid dynamics El ED EH El Ea. For surfaces, which will be our main focus, several authors have developed 
approaches to find approximate conformal parametrizations with respect to the unit disk or sphere |33l[3T||36l|26l[37j 
[68l[34l|27l[35l. More recently, quasi-conformal parametrizations based on the minimization of the Beltrami coefficient 
have been designed EHEDEtI- Another class of non-rigid registration methods include those based on optimal mass 
transportation (301 [32l [39l [40l, while (441 |42l [43l introduce comparison methods based on Gromov-Hausdorff or 
Gromov-Wasserstein distances. Computational methods based on integer programming and graph optimization have 
also been recently introduced (76l|78l[23l. We also refer the reader to the survey papers cmBiigHEa and textbooks 
(in[Til for additional entries on the literature. 

In this paper, we discuss an extension of the LDDMM framework, in which multiple shapes are registered simulta¬ 
neously within a deformation scheme involving contact constraints among the shapes. This is represented and solved 
as a constrained optimal control problem, in the spirit of the general framework recently introduced in (21 . 

Indeed, one of the characteristics of LDDMM is that it derives shape deformation from a global diffeomorphisms 
of the whole ambient space considered as a homogeneous medium, and does not allow for a differentiation of the 
deformation properties assumed by the shapes, or, more precisely, the objects they represent. This crude modeling may 
provide results that are not realistic in some applications. Consider the situation in which one studies several shapes. 
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representing, for example, different sub-structures of the brain. In this case, if one assumes that all shapes are deformed 
by a single flow of diffeomorphisms, shapes coming too close to one another will undergo a tremendous deformation, 
which creates artifacts that can mislead subsequent analyses. One would rather associate a different diffeomorphism 
to each shape, independent from the others, but the issue is that the resulting collection of diffeomorphisms may not 
be consistent: the shapes could overlap along the deformation. The solution briefly introduced in (Jl and developed 
in this paper is the following: embed the shapes into a ’’background”, complement of the shapes, deformed by a 
new, independent deformation, and add constraints such that, as all the shapes are simultaneously transformed, their 
boundary moves with the boundary of the background so that the configuration consistency is preserved. This is the 
approach that we develop here, focusing on surface registration. Note that a multi-diffeomorphism approach has been 
recently developed for image matching da, each diffeomorphism being restricted to a fixed region of the plane. The 
main (and fundamental) contrast with what we develop here is that, in our case, these subregions are variable and 
optimized, while they were fixed in da. The models along which sliding constraints are addressed in this paper and 
ours also differ. 

This paper is organized as follows. We start by recalling the classical LDDMM algorithm in Section setting 
the definitions, notation and appropriate framework for the rest of the paper. Then, in Section we introduce rig¬ 
orously the concept of multishape, describe identity and sliding background constraints, and describe the augmented 
Lagrangian algorithm for general constraints that will be used for in our numerical simulations. Section [^follows, 
specializing the algorithm to the case of identity and sliding constraints in great details. Finally, Section [^applies our 
method to synthetic examples to to real data as well. 


2. Large Deformation Diffeomorphic Metric Mapping 

2.1. Notation. In this paper, we define a shape as a embedding g : M ^ where M is a compact manifold. 
We denote by M the corresponding shape space, which is an open subset of the Banach space Q = C^{M, M^). 

Typical examples are as follows: 

• M = {1,..., m} is finite, and q can be identified with a collection gi,..., of distinct points in 

• M = [0,1] and g is a curve in 

• M = (the unit sphere in M^) and g is a hypersurface. 

Our goal is to discuss models in which several shapes can deform, while being subject to contact constraints. 
The deformation process will be similar to the one designed for large deformation diffeomorphic metric mapping 
(LDDMM), which can be formulated as an optimal control problem. Before introducing our general framework, it 
will be easier to start with a description of the now well explored single-shape problem upon which we will build. 
For this, we let II • Ik) be a Hilbert space of vector fields on assumed to be continuously imbedded in the 
space B := Co(M^,M^), the completion of the space of smooth compactly supported vector fields for the norm 
II • ||p,oo, which denotes the sum of supremum norms of derivatives of order p or less, with p > 1. Then V possesses a 
reproducing kernel, that is a mapping K : (x, y) ^ K{x, y), defined over x with values in the space of d x d 
matrices, such that all partial derivatives with order less than p with respect to each variable exist and 

K{-,y)aGV with {K{-,y)a, w)y = a^w{y), 

for all (a, y) G (M^)^. The LDDMM algorithm uses flows of time-dependent vector fields v{-) G 1/^(0,1; V). 

2.2. Registering Two Shapes Using LDDMM. 

2.2.1. General Problem. The general LDDMM problem is formulated as the infinite-dimensional optimal control 
problem consisting of minimizing the cost functional 

F{v) = \j\\u{t)\\ldt + U{q{l)\ 


( 1 ) 
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subject to the constraint 

dtQ(t) = u{t) o q(t) for a.e. t G [0,1], 

^( 0 ) = Qinit- 

This differential constraint is a control system, where the control is the time-dependent vector field ^(.) GL2(0,l;y), 
the solution of which is g(t, •) = Qiniti')) where ip is the fiow of diffeomorphisms generated by defined as the 

unique solution of the Cauchy problem dtip{t) = u{t) o ip{t), (/^(O) = id^^d. For every time t, we have ip{t, •) G Diff^, 
the set of p-times differentiable diffeomorphisms in 

The function f/ is a matching cost function, that is, a penalization that pushes the solution of 0-0 towards a 
target. It will be assumed to be Frechet differentiable from Q to M. To simplify the discussion, and because this covers 
most of the interesting cases in practice, we will assume that there exists some fixed measure um on M such that its 
derivative, denoted dU{q) or dUq when evaluated at g G Q, can be expressed in the form dUq = ZqUM for some 
(um- measurable) Zq : M meaning that 

\/heCP{M,R‘^), {dUg\h)= [ h{m)-Zg{x)dvM{x). 

Jm 

Throughout the paper, for any Banach space X, the notation {p\v) will be used to designate the application p{v) of 3. 
linear form p G X* to a vector v e X. 

Under these assumptions, one can prove that the gradient of the objective function F defined by O (which is a 
mapping on T^([0,1], V)) is given by 

VvF{u){t,^) - [ K{i,q{t,x))a{t,x)dvM{x), 

JM 

where K is the reproducing kernel of U, and a \ [0,1] x M ^ is a time-dependent function defined by a(l, •) = 
-Zqi^i) and 

(3) dtOL =—[du o qY'a. 

This result implies, in particular, that the solutions of 0-0 must satisfy the Pontryagin maximum principle (see 
ElllSSl), which is the following first-order necessary condition for optimality. Let Hu be the Hamiltonian defined by 

Hu{p,q) = {p\uoq)- i||«||U 

for every rt G U, every q £ Q and every p G Q*. If rt(') is an optimal control, solution of the optimal control problem 
then it must be such that 

(4) u{t) = argmax^i^^(p(t),g(t)), 
where (p, q) are solutions of 

idtq{t) = dpHu{p{t),q{t)), 

\dtp{t) = -dqHu{p{t),q{t)), 

(p is the so-called co-state, or adjoint state) and p(l) = —dUqpiy Indeed, it suffices to take p(-) = a^OM and to 
use properties of the reproducing kernel to check that all the conditions are satisfied. Moreover, Q then implies that 
u = q{x))a{x) doM at every time. 


2.2.2. Examples of shapes and matching cost functions. 
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Example 1. To start with a simple example, let M = {1, ..., m} so that a shape q = (g'(l), ..., q{m)) is a collection 
of landmarks, and consider the landmark matching cost function defined by U{q) = \q{^) ~ yk\‘^^ fixed 

y = iyi,...,ym) e in which | • | is the Euclidean norm on We then have 

m 

{dU,\h) = 2Y,{q{k)-ykfh{k), 

k=l 

or, to interpret this result in the general form provided above, dUq = zum with z{k) = 2{q{k) — yu) and um the 
counting measure on {1 ,..., m}. 

Example 2. If /i is a scalar measure on and z a /i-integrable M^-valued function defined on support (/i), we will 
denote hy zy the vector measure such that 


{zji \w) = I z ’ w dyU, 

J support(/x) 


where z • w denotes the standard euclidean dot product. 

Vector measures of the form zy are continuous linear forms over any space that is continuously imbedded in 
Co(M^,M^), and in particular over any reproducing kernel Hilbert space W. For such a space, equipped with a 
reproducing kernel x, the operator norm of z/i is given by 



zix) ■ {x{x, y)z{y)) dii{x) dij{y), 


and more generally, the norm of the difference between two such measures is 


\\z^i - zjl\\l = jj z{x) ■ {x{x, y)z{y)) dfi{x) dy,{y) 

^^^k{xix,y)ziy))dii{x)djl{y) + JJ z{x) ■ {xix,y)ziy)) dfl{x) djl{y). 


Note that W and V have no relationship one to each other, except that both have a continuous inclusion in Cq (M^, M^), 
so that X is different from K. 

One can deduce from this the surface-matching cost function introduced in ( 6 D, in which an oriented surface S 
is represented as a geometric current and a dual-RKHS norm between currents is used. Identifying surface currents 
with vector measures, this leads to the representation of S given by ys = where Ns is the unit normal to S 

and as its volume form. Assume that M (the parameter space) is an oriented 2D manifold so that S = q{M) is a 
surface, and let be a positively oriented volume form on M. For m G M let Nq{x) G denote the “area-weighted 
normal” to S = q{M) at q{x), defined by Nq{x) = dq^ei x dqxe 2 where (ei, 62 ) is an arbitrary basis of T^M such 
that 7 ^a;(ei, 62 ) = 1. Then 

(/is I u;) = / {woq- Nq) dy, 

Jm 

for every w e Cq (M^ , ) • 

Now, given a reproducing kernel x and a target surface S = q(M), we define the surface-matching cost by 

= \\yq(M) - yq(M)\\x- 
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Example 3. This cost function is actually a special case of the most general framework in which one compares compact 
/c-dimensional oriented submanifolds of which we briefly discuss hereafter. Given such a manifold, S, with a 
global parametrisation q : M ^ S, one can associate to any uj e Cq (M^, (A^M^)*) (the set of differential /c-forms 
on that vanish at inflnity), its integral 

(C5 |cj) = [ uj = [ q^uj, 

Js J M 


where denotes the pull-back of u on M. An RKHS, W, of such forms, is a Hilbert space continuously embedded 
in Cq (M^, (A^M^)*), with kernel x(x, y) taking values in the space of bilinear functions on A^M^ x A^M^. The linear 
form Cs then belongs in W*, and is a special form of a geometric current, as deflned in 1201 . If S = q{M) and 5 is a 
target manifold, they can be compared using the operator norm 

(5) U{q) = \\Cg(^M)-Cs\\%,. 

Now, if we consider q : M = (—1,1) x M^M^a smooth perturbation of q such that q^ = q-\- eSq -1- o{s) where 
qe{x) = q{s, x) for 5 G (—1,1) and x G M, we have, for = { 5 } x M, 


{dU,\5q) 






where £d/ds is the Lie derivative along the vector held ^ on M (which is equal to (1,0) G M x T^M at any location 
(e, x) G M) and 


( 6 ) 




with the isometry from W* to W. We next show that 

{dUq\6q)= aq’dqvolu-f ^q • dqvoldu, 

JM JdM 

where voIm and voI^m are the positive Riemannian volume forms on M and dM, and : M ^ (resp. Pq : 
dM M^) is such that aq{x) is normal to S = q{M) (resp. I3q{x) is normal to dS = q{dM)) at q{x). Using the 
Cartan magic formula we get 

^djdeiQ '^djdsd^q d{iQ j q^{(i cc)) , 

so that, applying the Stokes theorem, 

{dUq I 6q^ / '^djdsiQ dec) + / '^ojOsiQ 

JMo JdMo 

where is/ds denotes the contraction operator. Note that, for ^ 1 ,..., G T^M, 

^a/a£(g*da;)(o,^)(6, • • • , 6) = dccg(^) {Sq{x), ..., dg^C/e) 

= dujqix){Sq'^{x), dq^^i,dq^£,k), 


where 5q^{x) denotes the projection of 5q{x) on Tq(^j.^S^ (since the form vanishes if 5q{x) G Tqf^^^S = Tq(^j.^q{M)). 
Since the set of /c-forms is a one-dimensional space on M, this means that we can write, for some function aq such 
that aq{x) G (Tg(a,)5')^, with x e M, 


^a/a£(^*dcu)(o,^) = {{Sq • aq) yoIm)x • 


Similarly, 

^d/ds{Q'^^q)iO,x) = {{dq • Pq)yo\QM)x^ 

forpqix) G {Tq^,)dS)^. 

The data term deflned in 0 is derived from this general construction, using the fact that two-forms in (or 
(d — l)-forms in M^) can be identifled with vector flelds via cca,(ei, 62 ) = w{x) • (ei x 62 ). The current Cs is then 
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identified with the vector measure /i^. The form uj introduced in ([^ becomes, introducing a global parametrizaion 
g : M ^ 5 of the target S, the vector field 

jRd 

With this identification, we have aq = div(u;) o qNq where Nq is the oriented “area-weighted” normal to S at q{x) 
defined previously, and Pq = Tq x w o q, where Tq is the oriented “length-weighted” tangent to q{dM) given as 
Tq{x) = dqxCi, where ei is the unit positively oriented tangent vector at x along dM. 

Example 4. Returning to surfaces, the discrete case, in which triangulated surfaces are compared, is, for practical 
purposes, even more important. We here also consider the case M = ,m}, with an additional family F of 

facets, which are ordered triples (i, j, k) with i, j, /c G M. (We assume that F is a consistent with a manifold structure: 
The set Vi of indices that share a facet with i must form a chain and no pair of indices can be included in more than 
two facets.) 

Given a one-to-one mapping q : M ^ define Sq as the collection of triangles Sq = {(g(i), q{j), q{k)), (i, j, k) G 
F}. If / = {i,j,k),l&tq{f) = iqii),q{j),qik)),Ng{f) = (g(j) - g(i)) x {q{k) - q{i)) andc,(/) = {q{t)+qij) + 
q{k))/3 respectively denote the triangle, area-weighted normal and center associated to the facet /. Following lIMl . 
we define the vector measure associated to q by 

yq = 

feF 

Here, (with x G denotes the atomic measure of mass 1 with support {x}. The (discrete) surface matching cost 
associated to a target q is then defined by 

Uiq) = Wiig - 1,41. 

Note that q does not need to be consistent with g, and can be defined on a different set of indices, M = m} 

and triangle structure F. One then has dUq = aqUM, where, as above, um is the counting measure on M, and 

ag(i) = E (dZl^fFgif)/3 + egif,i) x Z(c,(/))), 
feF,ief 

with eq{f,i) = q{k)—q{j) the oriented edge opposed to g(i) in g(/), and 

Z{-) = E x{;Cg{f))Ng{f) - E x{;c4f)N4f). 

2.2.3. Reduced Problem. Since the optimal control must satisfy 

(7) v= K{-,q{t,x))a{t,x)duM 

Jm 

for some function a defined on M, it is natural to parametrize vhy a and use this function as a new control. We define 
the inner product 

{a , P) = / a{x)^K{q{x), q{x))P{x) duM{x) duM{x) 

Jmxm 

between two measurable functions a and P defined on M. If v is given by 0, the reproducing property of the kernel 
implies that bfv = ll«ll?. The optimal control problem ([T])-([^ is then equivalent to the reduced problem consisting 
of minimizing the cost functional 


( 8 ) 
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subject to the constraint (control system) 

(9) dtq{t,x)= / K{q{t,x),q{t,x))a{t,x)diyM{x), 

JM 

almost everywhere over the time interval [0,1]. 

According to EiiniisiiiTa, we have \/F{a) = a — p, where p (the co-state) is a time-dependent vector-valued 
measurable function on M such thatp(l)uM = — and 

dtp{t) = -dq (^{p{t), a{t))^ - ||aW||^(t)/2) , 

where q is defined by Here, the gradient is computed with respect to the inner product (•, 

3. Multiple Shape problems 

3.1. Motivating Examples. In the previous formulation, the shape evolution was controlled by a single, smooth 
vector field v, inducing a single diffeomorphism of restricted to the considered shape. This approach has been 
successfully used to model variations of single, homogeneous shapes, and led to important applications in compu¬ 
tational anatomy, including, among many other examples, the impact of pathologies like Huntington disease ca, 
schizophrenia l53l . and Alzheimer’s disease iMiiiaiMi on brain structures. This deformation model, however, is not 
well adapted in situations in which several shapes interact, or situations in which shapes have heterogeneous parts. Let 
us review some motivating examples. 

(1) Consider a schematic representation of a kite, or a manta ray, composed with a two-dimensional surface, 
representing the body, and an open curve attached to it representing the tail. When comparing two such 
objects, the body is assumed to only show small differences in shape, while the tail can vary widely. 

(2) Consider a two-dimensional representation of a mouth, with two curves representing the upper and lower lip. 
Because the mouth can be wide open or closed, it is not possible to consider its deformations as resulting from 
the restriction of a smooth diffeomorphism of M ?. 

(3) Finally, it is natural, when analyzing multiple organs in the human body, to consider multiple shapes, each of 
them being relatively stable (only subject to small deformations) while their position with respect to each other 
is subject to larger variations, so that the background (the intersection of their complements) is subject to very 
large deformations. Here again, modeling the whole process with a single diffeomorphism is not adequate. 

These examples suggest using multiple deformations applied to each component of the considered model. Gener¬ 
alizing 0- Q, consider parameter spaces Mi,..., for an n-component model. Each shape, or component, is a 
mapping q^^ G Qk • Mk The shape space will then be Q = Qi x • • • x To each shape, associate a 

control Uk ^ Vk, where Vk is an RKHS embedded in Cq (M^, M^) with the state evolution equation dtq^^^ = Uk o q^^\ 
We can then choose each Vk according to how “wildly” we want to allow the k-th shape to deform. These evolutions, 
however, must be consistent with each other, implying contact constraints that we will consider in two forms: 

• Identity constraints: These are constraints that make a subset of the k-th shape stay stitched to a subset of 

the l-th shape, so that these subsets coincide in and move identically along the deformation. Given some 
pair {k, 1) G n}^, and given a one-to-one mapping pki : Aki C Mk ^ Aik = gki{Aki) C Mi, one 

has q^^\x) = q^^\gki{x)) for every x G Aki. 

• Sliding constraints: These are constraints that force a closed submanifold of the k-th shape to slide on a 

corresponding submanifold of the l-th shape along the deformation, for all (/c, /). Here, we assume that all 
parameter spaces are orientable differential manifolds, and that all q^^^'s are immersions. Given some pair 
(k^l) G and a closed submanifold without boundary Aki G Mk, there exists a diffeomorphism 

Qki : Aki Aik G Ml onto a fixed closed submanifold Aik of Mi such that (x) = q^^^ {gki{x)) for every 
X G Aki- 
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Let US turn back to the examples mentioned at the beginning of the section. For Example (1), we can take Mi = 5'^ 
and M 2 = [0,1], and, letting xq represent the north pole in 5'^, impose We can then assign 

different deformation models to and via the metrics on Vi and V 2 . 

Example (2) requires a slightly more complex construction, that only imperfectly addresses the issue. Let Mi = 
M 2 = [0,1] and M 3 = {1, 2} x [0,1]. Let g^^^ represent the upper lip, g^^^ the lower one and g^^^ their union. We use 
the identity constraints g^^^ ( 1 ) = g^^^ ( 0 ) and g*^^^ ( 1 ) = g^^^ ( 0 ) for the extremities of each lip, and q^^^ ( 1 , •) = g*^^^ (•), 
g(3)(2^.) = g^^^(-). We take Vi = V 2 and choose V 3 such that the latter allows for large deformations at a small 
cost. With this model, it becomes easier to “almost” close the mouth, although the deformation inside the mouth must 
remain diffeomorphic, so that the closing cannot go all the way. 

For Example (3), there are n shapes, n — 1 of which are associated with the organs, and the last of which represents 
the background. For example, we can take Mk = for k = 1,..., n — 1, and Mn = {1,..., n — 1} x S‘^. 
Assuming that the shapes do not intersect, we can define identity or sliding constraints for the background, enforcing 
q^'^^ {k, •) = q^^^ (•) or q^'^^ {{k} x 5'^) = g*^^^ (5'^) = Sk for k e { 1 ,..., n — 1 } during the deformation. 

3.2. Induced Constraints. The previous constraints can be reformulated as equality constraints involving the state 
and control. Identity constraints q^^\x) = q^^\g{x)) are equivalent (taking time derivatives) to Uk{t^q^^\t^x)) = 
ui(t^ q^^^ (t, g{x))) as soon as the constraints are satisfied at time t = 0 , which we obviously assume. 

Making the same assumption, sliding constraints can be expressed as 

(10) N^^\t,q^^\t,x))'^{uk{t,q^^\t,x)) - Ui{t,q^'"\t,x))) = 0, 

where g^^^) is a d x [d — dim(A/i;/)) matrix consisting of independent vectors perpendicular to Tq(k)Bki{t) 

(e.g, a normal frame to Bki), with Bki = g*^^^ {Aki) for every (/c, /). Let us briefly justify this statement. 

We express the sliding constraint as q^^\t^x) = g{t^x)) for some diffeomorphism g{t, •) : A^i Aik, 

assuming a differentiable dependency on time. Taking time derivatives, we get 

Uk{t,q^'"\t,x)) = ui{t,q^’-\t,g{t,x))) + dq^’‘'> {t, g{t, x))dtg{t, x), x G 
Since g^^^ (t, g(t, x)) = g^^^ (f, x), we obtain 

Uk{t,q^'"\t,x))-ui{t,q^'"\t,x)) = dq^'‘\t,g{t,x))dtg{t,x) 

= dq^'^\t,x)dg{t,x)~'^dtg{t,x), 

which is tangent to Bki at x. Note that, since the image of •) is Aik for every time t, we do have dtg{t,x) G 
Tg{t,x)Aik = dg(t,x){TxAki), so dg{t,x)~^dtg{t,x) is well-defined. 

Conversely, assume that ( fT 0 | holds for every x G Aki, with q^^\0,x) = q^^\0, Qoix)) for some diffeomorphism 
go : Aki Aik C Ml. Then for every time t, the mapping 

w : X e Aki^ dq^^\t,x)~^{uk{t,q^^\t,x)) - Ui{t,q^^\t,x))) G T^Aki 

" -V-^ 

defines a time-dependent vector field on Aki • Since Aki is a closed manifold, this vector field is complete, and we 
denote its flow by /i(t, •) : Aki ^ki- Then 

dtq^^^ (t, h{t, x)) = Uk (t, q^^^ (t, h{t, x))) + dq^^^ (t, h{t, x))dth{t, x) 

= ui{q^^\t,h{t,x))), 

where the last identity holds for every x G Aku so that g^^^ (t, h{t^ x)) and g^^^ (t, go{x)) satisfy the same differential 
equation with the same initial condition and therefore coincide. Hence, letting g{t,x) = go{h~^{t,x)), we obtain 
q^^\t,g{t,x)) = q^^\t,x) for every x G Aki. 

It is possible to extend this construction to the case where Aki is a compact manifold with boundary. In this case, 
the matrix (t) must consist of a normal frame along dBki (t) and possesses therefore an extra column, and we get 
an additional constraint along the boundary of Aki • 
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We will need the constraints to depend smoothly on q, and therefore we will need a smooth representation of the 
normal space to q (a smooth map q^ N{q))m order to be able to use When this is not possible (or convenient), 

one can also use the alternative approach of introducing a new state, say (f, x), evolving according to 

( 11 ) = -duk{q^^^f 

which ensures that (t, x) remains perpendicular to as soon as this holds true at f = 0. The constraint 

x)^{uk{t^ x)) — ui{t^ x))) is now a smooth function of the extended state. 

The two problems that we consider are therefore special cases of the general problem considered in da), which is 
the problem of minimizing the cost functional 

n n 

( 12 ) 2 ^/ \Wk{t)\\ljt + Y,Uu{q^^'>{l)), 

^ k=l 

subject to the constraints 

(13) dtq^'^\t) = Uk{t,q^''\t)), and C{q{t))u{t) = 0, 


almost everywhere over the time interval [0,1], where C : M ^ L{V,y) takes values in the space of bounded linear 
operators from y to a Banach space y. Here, we have V = Vi x - • x Vn and q = {q ^^^,..., 

The study of this constrained optimal control problem, and in particular, the derivation of its first-order optimality 
conditions (of the type of Pontryagin maximum principle), is challenging in this infinite-dimensional setting. In (21, 
it is proved that, under some differentiability conditions, and under the important assumption that C{q) is surjective 
for every q ^ A4, optimal solutions must be such that there exist p = {p^^\ ... G H^{[0,1],Q*) and A G 

L2([0,1],3^*) that satisfy 


(14) 


'dtq^^^ =Uk{q^^^), 

= -dg(k) I - dg(k)(X I C{q)u), 

^ {uk,v)vk=-(p‘'^'’\voq^^'^'^-{X\Ck{q)v), vGVk, 

n 

'y2^k{q)uk = 0 , 

^ k=l 


where Ck{q)uk = C{q){0 ,..., 0,14/., 0 ,..., 0). 

Unfortunately, the constraints C{q) that correspond to our identity or contact constraints are, in general, not sur¬ 
jective, and the results of O cannot be applied in a fully general infinite-dimensional context. However, surjectivity 
becomes almost straightforward when these constraints are discretized to a finite number. They are true as soon as the 
points involved in the constraints are all distinct, which is a mild assumption. We now proceed to the description of a 
discrete version of this approach. 


4. Discrete Approximations 

4.1. Augmented Lagrangian. As an example, and to simplify the presentation, we detail our implementation for 
multi-shape problems in which shapes interact (through constraints) with a background, but not directly with each 
other. Direct interactions between shapes can be handled in a similar way. Our constrained optimization method uses 
the augmented Lagrangian method (see, e.g., 1501 ). In a nutshell, in order to minimize a function u F{u) subject 
to multi-dimensional equality constraints C{u) = 0, the augmented Lagrangian method consists of considering the 
functional 

Liu)=F{u)-X-Ciu)+'^\C{u)f, 
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in which A lives in the dual space of the space of constraints Y, and /i is a positive real number. Each iteration of the 
algorithm consists in minimizing L with fixed A and fi (our implementation using nonlinear conjugate gradient) until 
the gradient norm passes below some upper bound, and then in updating A according to the rule 

A ■{— A — 

before running a new minimization of L. The constant /i is increased only if needed, i.e., if the norm of the constraint 
did not decrease enough during the minimization. More details can be found in ISOl . 

We first apply this to identity constraints, which only require the shapes to be discretized into a sets of points. 
We will then discuss sliding constraints, which will require more structure in order to define normal frames to the 
boundary. 

4.2. Identity Constraints. We consider n — 1 objects, discretized into point sets, so that M/. is a finite set of indices 
for each k. Let and = {x^i \ ...,for k = 1,...,n — 1, with rrik = \Mk\. We add as 

n-th object the background, defined on = ({1} x Mi) U • • • U {{n — 1} x M^-i). We let = q^'^\k,j), 

= {z^\ ..., Zml) and ^ = {z^^\ ..., (a collection of m = = mi + ... + rrin-i points). 

Assume that end-point cost functions ), • • •, ) are defined, typically measuring the discrepancy 

between each collection of points and an associated target. We assume similar functions Ui{z ^^^),..., 
for the background, typically using Uj = Uj. The associated constrained optimal control problem consists in mini¬ 
mizing the cost functional 


-i n n—1 n—1 

k=l k=l k=l 

subject to the constraints (almost everywhere along [0,1]) 

'dtxf^ = Uk{xf^) j = l,...,mk, k = l,...,n-l, 

< dtzf'' = Un{zf'’), 3 = 1,.. .,mk, k = l,...,n-l, 

^z^^^ = x^^\ /c = l,...,n — 1. 

For y and y' ordered families of points in let {y^y') be the matrix formed with all d x d blocks Ky,, {yi ^y'j), 
Sindlot K^^\y) = K^^\y,y)Jork = 1 ,... , n, where is the kernel of 14. Since the problems only depend on 
the values taken by i^i,..., on their corresponding point set trajectories x^^\ ..., x^'^~^\z, the optimal vector 
fields take the form 

Uk{-) = K«(-, = 1,. .. ,n - 1, 

uy) = K^^\;z)p, 

for some families of d-dimensional vectors ..., The problem can therefore be reduced to the finite¬ 

dimensional optimal control problem consisting in minimizing the cost functional 

1 1 

E{a,l3,x,z) =-Y, - p ■ {K^'^\z)l3) dt 

2 Jo 2 Jo 

n—1 n—1 

k=l k=l 
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subject to the constraints (almost everywhere along [0,1]) 

dtZ = K^^\z)p, 

^{k) _ ^{k) ^ /c = 1 ,..., n — 1. 

Extending E with the augmented Lagrangian method, we introduce coefficients k = 1,..., n — 1 (where 
has the same dimension as x^^^) and /i > 0, defining 


n-l .1 n-1 

L{a,l3,x,z) = -V] / (A'')Vj7fe(:c('=)(l)) 

^ fe=i -^0 ^ -^0 fe=l 

n—1 n—1 n—1 

fe=l fe=l -^0 ^ fc=l -^0 

which will be minimized subject to the constraints (almost everywhere along [0,1]) 

\dtZ = K^^\z)p. 


From the constraints, L can be considered as a function of a and /3 only, and its differential with respect to these 
variables can be computed via the adjoint method as follows. Denoting the co-states by /c = l,...,n — 1, and 
, the associated Hamiltonian is 


n-l .1 .1 

k=i -^0 


The computation of the gradient of L follows the same general scheme as the one described in Section |T2] for the basic 
LDDMM algorithm. Given a and P and the associated trajectories x and 2 ;, one has solve the adjoint equations 

dtP^ = -d,H, p"A(i) = _vi/fc(^W(i)), fc = l,...,n-l. 

The computation of the differential system gives 


• • ^ K ' ‘ ^ K 

dtPT" = -E^iipT'' ■ K^'^\xf\xP)aP) - y] Vi(af) • K^^\xf\xf)pf) 


i=i 


i=i 


E2Y,^i{oii ' - (A^^) 

i=i 


and 

dtpf = 


n —1 mi 


n—1 mi 


1=1 j=i 1=1 j=i 

n—1 mi n 

1=1 j=i 1=1 
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The gradient of L with respect to {a, /3) is then deduced from the partial differentials of H with respect to these 
variables, yielding 

= xW(A''))(«('=) 

VpL = K^^\z){/3-p^). 

Alternatively, on may choose to use the gradient relative to the dot product on x ... x I 4 , which is simply given 
by 

VpL = /3-p^. 

The latter choice is simpler, and generally more efficient numerically. 


4.3. Sliding Interface. Assume that the parameter sets Mk are vertices of pure oriented, simplicial complexes Tk of 
dimension < d (we will however only provide implementation details for codimension d — = 1). We let 

denote the set of facets of the k-th complex. We also assume that Ti,..., T^_i are disjoint and that Tn is their union, 
Tn = Ufc=i We also let F = [JlZ\ Fk (disjoint union). 

The associated shape space is formed by functions : Mj^ such that Qkif) is not degenerate (i.e., has 

maximal dimension) for all f e Ff^. Each object is allowed to slide against the background. We will write = 
q{k) (^Mk), /c = 1 ,..., n — 1 , and {Mk), z = {zM ^^ in accordance with our previous notation. 

If / G F/c is a facet in Tk C Tn, we discretize into 


(15) 


7V(")(/). 



I = 0, 


where N^^\f) is a d x {d — Vk) matrix spanning the normal space to q^^\f), assumed to be defined as a smooth 
function of q^'^\ If Vk = d — 1, this is always possible, since is a vector that can be taken as the cross product 
of Zf ^2 — • • • 5 — ^f,i where ..., Zf^d is any labeling of the vertices of q^'^^ (/) ordered consistently with 

the orientation. 


We now restrict to this case, with d = 3, so that shapes are triangulated surfaces in as discussed in Section 2.2 
For f G Fk and j G /, we denote by ejj the edge {zj^) — zjf^) where j' and j" are the other two vertices of / such 
that (j,is positively oriented. Similarly, let e' j = (zjf^ — Zj^^) and = {zj^} — z^^^) be the two edges 

(k) 

Stemming from Zj so that 

X e"^ = 2area(</(")(/))Ar(")(/) =: N^^\f) 

is the area-weighted positively oriented normal to / in q^^^ {Mk). Note that ejj = e'jj — e’-j. 

With this notation, we can rewrite the constraint in the form 


y]det(e'J, e'-j,Uk{zf^) - Un{zf'^)) = 0. 


holding for all f G Fk and /c = 1 ,..., n — 1 . 
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Introducing a Lagrange multiplier A/ for each of these constraints, after reduction of the vector fields, which 
proceeds similarly to the identity constraints case, the augmented Lagrangian takes the form 


L{a,P,x,z) 


with 


^ k=l k=l 

n—1 n—1 

k=l k=l feFk 


-p{k) 


{x^^\z) ;= ^det {e'^j,elf,K^^\zf\x^^">)a^'^'^-K^'^\i 


ik) 



We now compute the evolution equations for the co-states, as done with identity constraints. For f e Fk and 
i G M/c, we have 

(16) ^Vi(af) ■ K('^\xf\zf^)N^^\f)). 

Denoting 

■= - Un{z^’'^)), 

j^f 

ifiefe Fk, then 

n—1 mi 

(17) d^wrf = -eijx6^'^\f)- EEv 1 (/) • ) 

1=1 3 = 1 

mk n—1 mi 

+ ^Vi(7V(")(/) )af) - Vi(/3P ■ K^^\zf\zf)N^^\f)). 

3=1 1=1 3=1 

Let _ _ _ ^px,n-i gjjjj pz _ (j)Z,i^ _ _ _ be the co-states. Let = A/ — ^ let 

Fkii) = {/ e Mfc : i e /}. Then 


dtPT'^ = 


Nk 

E 

i=i 


Nk 


- • K(’"\xf\xf)af'^) - 51 Vi(af) • 

i=i 

Nk 

+ 2 5^Vi(afLifW(a;W,a;W)«('=))- 5] 7Wa^(.)rf, 


j = l 


/GFfe 


and 


AT 


AT 




i=i 


i=i 


AT 


+ 25^ Vi(/3P ■ K^^\zf\z^)(i^) - Y, 

i=l /GFfc 


where d^{k)T^p and are given by ( p^ and ( [Tt] ). 
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For / G F/c (/c = 1,..., n — 1), we have 

j^f 

jef 

Letting 9j’^^ = (/)> the gradient of L in a and in fj is then given by 

V^(k)L = (A''))(«('=) - K^^\x^'‘\z)e, 

VpL = K^^\z){/3 - p^) + K^'^\z, z)e 

or, taking the Hilbert gradient, 

V„wL = 

VftL^ p-p^ + K^'^\z,z)e. 

In spite of it requiring the inversion of a linear system in the first equation, we found the latter version preferable to 
the gradient in our experiments. 

4.4. Remarks. 


Existence of constrained solutions. It is important to note that, according to in Theorem 1], there always exists at 
least one solution of (p^-(p^ satisfying the constraints. 


Convergence to surfaces. A question naturally arising is whether our discrete approximation using triangulations 
converges to the smooth setting as triangles get smaller and smaller. More precisely, assume that smooth initial 
surfaces are triangulated, with increasingly fine triangulations = 1 , 2 ,..., 

where labels the vertices of a simplicial complex whose faces are We discuss whether minimizers 
{uk/^ /c = 1 ,..., n) of the discrete problems have a subsequence that converges to a minimizer {uk^ /c = 1 ,..., n) of 
the limit problem. 

Assume that the following condition holds for the sequence of triangulations: 


(i) We assume that for all k and i, and for every / G Fk^i, there exists an embedding ^ ^ 

fk) 


-f-f V q{k) _ 
^k,i ^ ^init — 

(where is the interior of the triangle quAf)) that / G F/^,^) partitions 

up to a negligible set and maxj ^ — id^^ ^ ||i,oo ^ 0 when I ^ oo. 

This conditions ensure that data attachment terms like those described in section |2.2.2| computed at diffeomorphic 
transformations o converge to the same term computed at cpk o as soon as converges to in 
Given this, we sketch the argument leading to the consistency of the discrete approximations. 


For identity constraints, one can use O Proposition 5], which proves that, if the triangulations are nested (every 
vertex at step i lies on the limit surface and is also a vertex at step ^ + 1), then one can extract, from a corresponding 
sequence of identity-constrained optimal vector fields, a subsequence that converges towards an identity-constrained 
solution of ([T^-([Tj|). 

For sliding constraints, one cannot directly apply O Proposition 5], because the constraints are not nested, even 
when the triangulations are. To obtain a consistent approximation, we need to relax the discrete problems. More 
precisely, let t i-G^ ... ,Un{t)) G Vi x • • • x Vn be a minimizer of the continuous problem with sliding constraints, 

and let (cpi,... ^cpn) denote the corresponding fiow with q^^\t) = (pkif) o q^^.^ the corresponding deformation of 
g-^]^, and ^ normal to Sk{t) = q^^\t){Mj^) at q^^\t,x). In particular we have, for 

every /c = 1,..., n, every x G Mk, and almost every time t, 

N^^\t,x) ■ {ukit,q^^\t,x)) - Un{t,q^'"\t,x)) = 0. 
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Moreover, as ll^/c(^)|P is constant, both Uk{t) and duk{t) are a-Lipschitz for some positive constant a that 

does not depend on t or /c. 

Now let = (pk{t) o be the corresponding deformation of the discretization at step i. Recall that 

Nik,e) denotes the unit normal to the triangle (t, /). We will prove: 

(ii) The discretized deformations at step £ satisfy the following relaxed sliding constraints 


EK(-i )-«n(.r')| <e. 

\jef 

at almost every t and for every face / in for a suitably chosen sequence ££ > 0 going to 0 as ^ goes to 


Indeed, fix a face / and an integer £. Define 

for every time t. Note that assumption (i) implies that 


for some sequence 7 ^, independent of / and t, and going to 0 as ^ goes to infinity. Assumption (i) also implies that 
there exists a sequence 7 ^, independent of / and going to 0 as ^ ^ oc such that for y G ^(/), and = y, we 

have 

/)| + E l9St(^) - 

j£f 

Some triangle inequalities, Gronwall’s lemma, and the fact that and du^ are a—Lipschitz then imply 





< e£, 


with ££ = ar]ie^ + 7 ^ going to 0 as goes to infinity. 

Consequently, if ( 1 ^ 17 ,... is a sequence of minimizers of the discretized problem at step £ with relaxed 

sliding constraints 




\jef 


j 


{t)) - Un,i{t,Z. 


ik,i) 



< Si, 


we see that the infimum limit over £ of the respective discretized costs of ..., is smaller than or equal to the 
cost of a minimizer of the continuous problem with sliding constraints. So to prove that a limit point of that sequence 
is a minimizer of the cost for the continuous problem with sliding constraints, all we need is to check that any such 
limit point does satisfy the constraints. 

So let ..., Un^i) be a sequence of minimizers of the relaxed discrete problem at step ^that weakly converges 
to (r^i,..., in 1^1 X • • • x Vn (which is true for at least one subsequence of any minimizing sequence), then the 
associated flows ipk.i and their first two space derivatives converge uniformly in time and space to the flows (pk and 
their first two derivatives. 

From this it is easy to see that any sequence approximating as in assumption (i) is such that (f, j^) 

q^^\t,x) and at all times. Moreover, such a minimizing sequence must be, like all 

geodesics, such that Ylk=i \Wk,e{t)\\y^ is constant in time, and smaller than the cost function associated to, say, 
Uk = 0 for all k. This implies that the vector fields Uk/{t) are continuous uniformly in k, £ and t, which, combined 
with the continuity of the evaluation functionals in an RKHS implies that ji)) Uk{q^^\t^ x)) for all 


fit is easy to prove that such minimizers exist using the same method as that of (2), and replacing equality constraints with inequality constraints. 
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times. Consequently, one easily checks that each , Un,i) satisfy a relaxed version of the continuous version 

of the constraints, with a precision that goes to 0 as ^ goes to infinity. This finishes proving that the constraints are 
satisfied with exactitude at the limit. 

Kernel derivatives. Expressions similar to Vi (n • if (x, y)a) appear at multiple times in the previous computation (for 
some vectors n and a). For radial kernels (K{x,y) = G{\x — ^p)Id]^d), we have 

Vi(n • K{x, y)a) = 2G'{\x - y\^){n ■ a){x - y), 

which (slightly) simplifies the expressions. 

Sliding Interface - Alternate Version. As discussed in Section the sliding constraint can also be handled by intro¬ 
ducing a new state variable N that tracks a vector (or frame) normal to the interface via O- In the discrete case, one 
can discretize this equation by introducing states N{f), f G F, indexed by the facets of M, and evolving according to 

ief 

where |/| is the number of vertices in /. The sliding constraints are now expressed in terms of the state variables 
in a more direct way, but with a new co-state variable for the normals, bringing in an extra degree of complexity 
and increasing the computational cost. Note that the finite-dimensional reduction is still possible in this case, so that 
^(n)(.) _ evolution of the normals can be expressed in a form involving the differential of the 

kernel. This yields an adjoint system involving second derivatives of the kernel. We will not detail the computations in 
this paper, since they follow the same pattern as the other two that were already discussed (see l(56l for more examples 
on how higher-order variables can be handled in similar contexts). Note that this alternate version of the sliding 
constraints is slightly more general than the one discussed in the previous section, since it does not require a definition 
of a normal field that smoothly depends on the manifolds. 

5. Experimental Results 

5.1. Synthetic Example. The first example is described in FigureIn this synthetic example, the template has two 
identical balls initially close to each other. In the target, the first ball (referred to as “Ball A”) gets bigger, and “impacts” 
the other one (referred to as “Ball B”), which assumes an oblong, non-convex shape (the target shapes slightly overlap, 
so that an exact homeomorphic match cannot be achieved). 



Figure 1. Template and target shapes for synthetic example 
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Our results, provided in Figurestoillustrate our multishape deformation method, and use two complementary 
deformation indexes: 

(i) the tangent Jacobian, which is the Jacobian determinant of the surface-to-surface transformations, and which 
measures the ratio between the areas of elementary surface patches at each point before and after deformation; 

(ii) the normal Jacobian, which is the ratio of the Jacobian determinant (of the 3D diffeomorphism) to the tangent 
Jacobian, and which measures the ratio between the length of an infinitesimal line element normal to the 
surface after and before transformation. 

These indexes are mapped on the deformed template image, which is close to the target. 

Figure [^compares the normal Jacobian of the shape and background deformations when using identity constraints. 
While shape diffeomorphisms characterize each shape transformation (uniform variation for Ball A, expansion at the 
top and compression otherwise for Ball B), the effect of compressing the space is clearly visible in the background 
deformation, when the two shapes get close to each other. 

Figure [^provides the corresponding tangent Jacobian, which is identical for shape and background transformation 
since we are using identity constraints. 

Figures 1^ andcompare the normal and tangent Jacobians for the synthetic experiment with sliding constraints. 
Regarding the former, the most notable difference is with Ball B, which shows an expansion pattern at the tips in its 
shape diffeomorphism which is inverse of the one observed with identity constraints. One plausible explanation is that 
sliding constraints allow the two shapes to use translation-like motion to position themselves differently, without the 
need for limiting the amount of shear in the background that would have resulted from identity constraints. The second 
notable difference can be noted in the background diffeomorphism, in which compression is mostly observed with Ball 
B. In contrast with the identity constraints, the tangent Jacobians are very different between shape and background 
diffeomorphisms. Note that Figureuses two different color scales for the left and right panels because of the strong 
difference between the ranges of the Jacobians in each case. The background deformation, in particular, has a huge 
tangent expansion around the impact location, which cannot be observed in the shape deformations. Note that both 
patterns in the sliding case are very different from the one that was observed in the identity case. 

For comparison purposes. Figure [^provides the result of the LDDMM algorithm using a single diffeomorphism. 
One observes a very strong compression effect for the normal Jacobian resulting in an expansion observed on the 
tangent Jacobian on Ball B, that was not observed in any of the constrained examples. The nice uniform expansion in 
Ball A that could be observed in the sliding constraint case is not observed either. 


5.2. Subcortical Structures. We now describe an example mapping a group of three subcortical structures: hip¬ 
pocampus, amygdala and entorhinal cortex (ERC). The template and target sets are represented in Figure |7] One can 
observe shape changes in each structure, combined with a significant displacement of the ERC relative to the other two 
structures when comparing template to target. Because the structures were segmented independently, there is some 
overlap between the target hippocampus and amygdala. 

Figures and [^provide the normal and tangent Jacobian obtained with identity constraints, while Figures [T^ and 
[^provide this information for sliding constraints. The two types of constraints provide similar deformation indices, 
especially for the normal Jacobians (Figures and [T^. Minor differences in the tangent Jacobian can be observed 
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Figure 2. Three views of the normal Jacobian with identity constraints: shape diffeomorphisms 
(left) and background diffeomorphism (right). 
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Figure 3. Tangential Jacobian: shape and background diffeomorphisms (identity constraints). 


(Figures]^ and [TT]). The deformation patterns associated to using a single diffeomorphism (Figure[^ are significantly 
different, though, exhibiting very strong compression, for example, where shapes are close to each other. 

6. Discussion 

The previous approach provides a solution, using constrained optimal control, of the important issue of dealing 
with multiple objects with varying deformation properties for registration. We have focused on surface matching, 
numerically dealing with constraints using an augmented Lagrangian method. Note that a similar approach was 
introduced for plane curves in (Jl . 

The formulation is quite general and can accommodate constraints in various forms, including the examples dis¬ 
cussed in SectionThe investigation of these additional applications will be the subject of future work. One of the 
limitations of the present implementation is the slow convergence of the augmented Lagrangian procedure, for which 
each minimization step is, in addition, high dimensional and computationally demanding. One possible alternative can 
be based on solving the optimality conditions ( p^ (which hold in the discrete case) by means of a numerical shooting 
method. This approach has, however, its own numerical challenges, because solving ( p^ requires the determination 
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Figure 4. Three views of the normal Jacobian with sliding constraints: shape diffeomorphisms 
(left) and background diffeomorphism (right). 
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Figure 5. Three views of the tangent Jacobian with sliding constraints: shape diffeomorphisms 
(left) and background diffeomorphism (right). 
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Figure 6 . Three views of the Normal (left) and tangential Jacobians (right) when using a single diffeomorphism. 


of A such that the last equation (constraint) is satisfied, and this leads to a possibly ill-posed problem for systems in 
large dimension (see O for additional details). 

We have illustrated our examples using deformation markers derived from the Jacobian determinant. This markers 
are routinely used in shape analysis studies and led to important conclusion in computational anatomy. When dealing 
with multiple shapes, however, figures [6| and \V^ show that, when using the classical LDDMM method with multiple 
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Figure 7. Template (blue) and target (red) shapes for subcortical structures. The hippocampus is 
the central shape, with the amygdala on its left and the ERC on its right. 


shapes, these markers becomes as much, if not more, influenced by interactions between the shapes as by the changes in 
the shapes themselves. For this reason, multi-shape computational anatomy studies have applied registration methods 
separately to each shape, without ensuring that the obtained diffeomorphisms are consistent with each other. This 
limitation is addressed in the present paper, in which we exhibit deformation markers that are meaningful in describing 
tangential and normal surface stretching, while being consistently associated to a global transformation of the space. 
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